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Abstract 



An approximate analytical approach to describe the 
stochastic motion of sound rays in deep ocean is 
developed. This is done for a realistic propaga- 
tion model with an internal wave induced pertur- 
bation imposed on the smooth background sound 
speed field. The chaotic ray dynamics is analyzed 
using the Hamiltonian formalism taken in terms of 
the action-angle canonical variables. It is shown that 
even at ranges of a few thousands km, the magni- 
tude of range variations of the action variable is still 
small enough to be used as a small parameter in the 
problem. A simple expression for the difference in 
travel times of perturbed and unperturbed rays and 
approximate analytical solutions to stochastic ray 
(Hamilton) equations are derived. These relations 
are applied to study range variations of the timefront 
(representing ray arrivals in the time-depth plane). 
Estimations characterizing the widening and bias of 
timefront segments in the presence of perturbations 
are obtained. Qualitative and quantitative explana- 
tions are given to surprising stability of early portions 
of timefronts observed in both numerical simulations 
and field experiments. 

1 Introduction 



We consider the ray dynamics in a deep ocean acous- 
tic waveguide with an internal wave induced per- 
turbation to the smooth (background) sound speed 
field. It is assumed that statistics of random internal 



waves are determined by the empirical Garrett-Munk 
spectrum Jl], ||, ||. Numerical simulations demon- 
strate that although this perturbation is weak, it 
gives rise to a rather strong ray chaos || [|. In the 
presence of internal waves ray trajectories exhibit ex- 
treme sensitivity to initial conditions and at ranges 
of a few thousand km parameters of a ray trajectory 
with the given starting depth and launch angle are 
practically unpredictable and can be described only 
statistically. On the other hand, even under condi- 
tions of ray chaos the arrival pattern retains some of 
its features observed in the unperturbed waveguide. 
In particular, both numerical simulations || ||] and 
field experiments || || show that early portions of 
arrival patterns formed by steep rays manifest sur- 
prising stability. Due to this property, the ray travel 
time is the main signal parameter in the underwater 
acoustics experiments from which inversions is per- 
formed to reconstruct ocean temperature field Q], |J. 

Our main objective is to develop an approximate 
analytical approach for description of the ray motion, 
including analysis of ray travel times, at very long 
ranges. We argue that this goal can be accomplished 
in the scope of a perturbation theory based on the 
Hamiltonian formalism in terms of the action-angle 
variables ^ ftO[ . Range variations of the refractive 
index in our propagation model are not adiabatic and 
we cannot use the adiabatic invariance of the action 
variable. Nevertheless, it turns out that, even at long 
ranges, the variance of the action variable is small 
enough to be considered as a small parameter in the 
problem. 

In the present paper we derive an approximate an- 
alytic relation for the difference in travel times of two 



1 



rays one of which propagates in the perturbed waveg- 
uide and another one in the unperturbed waveguide. 
We also obtain approximate solutions to the stochas- 
tic Hamilton (ray) equations for the action and angle 
variables. These solutions are then applied to analyze 
statistics of ray travel times at ranges up to 3000 km. 

Our attention is focused on the so-called timefront 
representing ray arrivals in the time-depth plane. We 
estimate the width of different branches (segments) 
of the timefront in the presence of perturbation and 
establish the criterion of nonoverlapping of neighbor- 
ing segments. It is also shown that the perturbation 
causes not only a widening (dispersion) of timefront 
segments but some regular bias of the segments as 
well. The estimation of this bias is also presented. 
Predictions made with our approximate analytical 
approach are verified by comparison to results of nu- 
merical simulations. 

Note, that in the present paper we consider only 
internal-wave-induced perturbations. An important 
issue of travel time biases due to mesoscale inhomo- 
geneities (see, e.g., Ref. |ll||) has not been broached 
here. 

The paper is organized as follows. The ray (Hamil- 
ton) equations in terms of the position-momentum 
variables, (p, z), as well as an environmental model 
on which we rely in this paper are presented in Sec. 
|2|. The emphasis in this section is on discussion of 
properties of timefronts obtained by numerical solu- 
tion of the ray equations. The action- angle canonical 
variables, (1,0), are introduced in Sec. [|. Section f| 
presents the derivation of our main relation for the 
difference in travel times, At, of perturbed and un- 
perturbed rays. In Sec. || we deduce approximate 
stochastic equations governing fluctuations of action 
and angle variables due to random inhomogeneities. 
Solving this equation yields statistical characteristics 
of / and 9. In Sec. || these characteristics are used to 
analyze relative magnitudes of different constituents 
of At and significantly simplify the expression for At 
by neglecting small terms. Our final expression for 
travel time variations is applied to investigation of 
the timefront structure. This is done in Sec. 0. Sec- 
tion U is concerned with generalization of the analyti- 
cal approach to a more realistic environmental model. 
Our results are summarized in the final section. In 



the Appendix we shortly discuss how the transforma- 
tion from (p, z) to (I, 6) variables can be performed 
numerically using a standard ray code. 

2 Timefronts in the presence of 
internal-wave-induced pertur- 
bation 

2.1 Ray dynamics in terms of the 
Hamiltonian formalism 

Consider wave propagation in a two-dimensional 
medium with the coordinates r (range) and z (depth). 
It is assumed that the z-axis is directed downward 
and the plane z = is the sea surface. The ray 
trajectory z(r) is determined by the sound speed 
field c(r, z) and can be found from Fermat's principle 
JlO| , |l2] , [l3| , |l4[ according to which the first variation 
of the functional 

ds 



c(r, z) 



= dr n(r, z(r)) 




(1) 



vanishes at the ray trajectory. Here n(r, z) = 

c r /c(r, z) is the refractive index, c r is a reference 

/ 2 \ 1/2 
sound speed, and ds — dr ( 1 + (dz/dr) J is the 

arc length. The functional S represents the so-called 
eikonal and it is related to the ray travel time, t, by 



t = S/c r 



(2) 



Formally considering Eq. ([j]) as an action function 
of some mechanical system with the r-variable play- 
ing the role of time, one can apply the standard rela- 
tions of classical mechanics [l^, This yields 
explicit expressions for the momentum, 



dz/dr 



y/1 + (dz/drf 



(3) 



2 



and the Hamiltonian, 



H = -^Jn 2 -p 2 . (4) 

Equations p = n sinx and H — — ncosx relate the 
momentum and the Hamiltonian to the ray grazing 
angle, \ |j§. 

Expression (||) for the eikonal can now be rewritten 

as 

S = J (pdz - Hdr). (5) 

Ray trajectories are governed by the Hamilton equa- 
tions |, | 

dr ~ dp ~ ff ~ v^^' 1 J 

dp 9i? ndn/dz , N 

— = = — . (7) 

dr dz ^/n 2 - p 2 

Equation C0) present the Hamiltonian formal- 
ism in terms of the momentum-position canonical 
variables. Later on (in Sec. ||) we shall introduce an- 
other pair of canonical variables, namely, the action- 
angle variables. 

2.2 Environmental model 

In what follows we shall consider a model of the sound 
speed field in the form 

c(r,z) = co(z) + Sc(r,z), (8) 

where cq(z) is a smooth (background) sound speed 
profile, and 5c(r, z) is a range-dependent perturba- 
tion. 

The unperturbed profile cq(z) used in our numer- 
ical simulation is typical for deep water acoustic 
waveguides. It is shown in the left panel of Fig. 1. 
The sound-channel axis, i.e. minimum of the sound 
speed profile, is located at a depth of 0.738 km. 

We consider an internal-wave-induced sound speed 
perturbation, Sc(r,z), with a zero mean (< 
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Figure 1: Left panel: Background sound speed pro- 
file. Right panel: Sound speed perturbation versus 
depth at three different ranges. 

6c(r,z) >— 0) and assume that the statistics of 
the internal wave field is described by the empirical 
Garrett-Munk spectrum |Q. A numerical technique 
for generation of such perturbations developed by J. 
Colosi and M. Brown has been applied. Equa- 
tions (1) and (19) from Ref. [|| have been used to 
generate a particular realization of Sc(r, z) which is 
used throughout this paper. It has been assumed 
that the buoyancy frequency profile v(z) is expo- 
nential, v(z) — vq exp(— z/B), and determined by 
two constants: a surface-extrapolated buoyancy fre- 
quency Mo = 27r/10 min -1 and a thermocline depth 
scale B = 1 km. We consider the internal wave field 
formed by 30 normal modes and assume its horizontal 
isotropy. Components of wave number vectors in the 
horizontal plane belong to the interval from 27r/100 
km -1 to 2ir/4 km -1 . An rms amplitude of the pertur- 
bation, (Sc) lms , scales in depth like exp(— ?>z/2B) and 
its surface-extrapolated value in our model is about 
0.5 m/s. Depth dependencies of 5c at three different 
ranges are shown in the right panel of Fig. 1. 
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Figure 2: Early portion of the timefront at the range 3000 km without (upper panel) and with (lower 
panel) internal waves present. Identifiers of rays forming some particular segments are indicated next to the 
corresponding segments. In the upper panel, arrivals with identifier +124 are depicted by a thick solid line. 
In the lower panel, arrivals with this identifier are marked by thick points. 
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Figure 3: Late portion of the timefront at the range 3000 km without (upper panel) and with (lower panel) 
internal waves present. In the upper panel, arrivals with identifier +160 are depicted by a thick solid line. 
In the lower panel, arrivals with this identifier are marked by thick points. 
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2.3 Numerical simulation of time- 
fronts 

Figures 2 and 3 show early and late portions of the 
timefront at 3000 km ranges, respectively, for rays es- 
caping a point source set at a depth of 0.78 km. The 
timefront in the unperturbed waveguide graphed in 
the upper panels of Figs. 2 and 3 has been com- 
puted using a fan of 16000 rays with starting mo- 
menta equally spaced within an interval correspond- 
ing to launch angles ±12°. The timefronts in the 
perturbed waveguide has been produced by tracing 
49000 rays with starting momenta covering the same 
interval. 

The timefront in the range-independent waveguide 
has the well-known accordion-like shape consisting of 
smooth segments (branches) ^ [l5j. Each segment 
is formed by points corresponding to arrivals of rays 
with the same identifier ±J, where J is the number of 
ray turning points and symbols + and — correspond 
to rays starting upward and downward, respectively. 
So, we can associate each segment with the identifier 
of rays forming this segment. Identifiers for some par- 
ticular segments in the unperturbed waveguide are 
indicated in the upper panels of Figs. 2 and 3. It 
is seen that the travel time grows with J. This is 
a typical situation for a deep water waveguide |[T^| : 
steep rays usually have greater cycle lengths (smaller 
J) and arrive earlier than flat ones. Segments corre- 
sponding to rays with launch angles of the same sign 
form a broken line. Two such lines shifted along the 
i-axis form the unperturbed timefront (see plots in 
the upper panels of Figs. 2 and 3). 

A very interesting and important feature of the 
perturbed timefronts is a remarkable stability of seg- 
ments formed by early arriving steep rays. This prop- 
erty of steep rays is well-known and it has been ob- 
served in both numerical simulations and field exper- 
imentsilU. 

Figure 4 presents the ray travel time, t, as a func- 
tion of Xs, magnitude of the launch angle. Here and 
in the remainder of the paper we present numerical 
results only for rays starting upward (for rays start- 
ing downward the results are absolutely the same). 
In the range-independent case the function t(xs) is 
smooth and monotonous (solid curve). In contrast, 
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Figure 4: Ray travel times at the range of 3000 km 
as a function of launch angle, without (solid line) and 
with (points) internal waves present. 

for the range-dependent case, we observe strong sen- 
sitivity of ray travel time to starting angles: points 
depicting travel times of perturbed rays are randomly 
scattered. Another manifestation of stochastic ray in- 
stability in the presence of internal waves is seen in 
Fig. 5 where the ray identifier at 3000 km range is 
shown as a function of the launch angle. Perturbed 
rays with close initial conditions may have quite dif- 
ferent number of cycles. Although steep rays look 
less chaotic compared to flat ones (see Fig. 4 and 
the upper panel in Fig. 5), a magnified view of the 
angular interval \s > 9° shown in the lower panel of 
Fig. 5 demonstrates that the dependence of J on Xs 
for steep rays in the perturbed waveguide is also split 
into a set of irregular stripes. 

The time spread, St, for rays with close launch an- 
gles in the perturbed waveguide can be estimated as 
a width of the area occupied by randomly scattered 
points in Fig. 4. For rays with Xs < 9° this yields 
St 1.5 s. However, if we select a group of arrivals 
corresponding to rays with the given identifier, we 
shall find that the time spread of these arrivals will 
be significantly less than St. It is seen in Figs. 2 
and 3 where two such groups of arrivals are shown by 
thick points. 

So, in spite of strong sensitivity of ray parameters 
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Figure 5: Points in the top figure show the ray iden- 
tifier at the range 3000 km as a function of launch 
angle, in the presence of internal waves. The bottom 
figure presents a magnified view of this dependence 
for steep rays taken directly from the top figure. The 
solid line in the top figure graphs the dependence 
J(Xs) m the unperturbed waveguide. 



fan. In the presence of perturbation such rays appear 
(see horizontal stripes with J = 112, 113, and 114 in 
both panels of Fig. 5 at Xs close to 12°) because 
perturbation leads to widening of the interval of ray 
grazing angles. As it has been indicated already, in 
our waveguide steep rays arrive earlier than flat ones. 
So, it is natural that in the presence of perturbation 
we observe rays arriving earlier than those from the 
unperturbed fan. However, it is surprising and abso- 
lutely unexpected that steep rays which appear due 
to scattering at random inhomogeneities, form quite 
regular segments. Moreover, these segments coincide 
with the unperturbed segments missed in the upper 
panel of Fig. 2. 

One of our objectives is to explain how chaotic be- 
havior of ray paths is compatible with stability of 
early portions of the timefront formed by steep rays 
and with unexpectedly small time spread of clusters 
of flat rays. Later on, we shall address these issues us- 
ing the Hamiltonian formalism in terms of the action- 
angle canonical variables 0, The action-angle 
variables are introduced in the next section. 



to initial conditions, numerical simulations demon- 
strate an unexpectedly small time spread for rays 
with the given identifier. Loosely, we can state that 
although the ray travel time exhibits a chaotic and 
unpredictable dependence on the launch angle, its 
dependence on the ray identifier is much more pre- 
dictable. This effect is most apparent for steep rays 
which form segments of the perturbed timefront al- 
most coinciding with the corresponding segments of 
the unperturbed timefront. However, this fact does 
not mean that steep rays with the same launch angle 
in the perturbed and unperturbed waveguides follows 
close ray paths. 

To illustrate this statement and demonstrate that 
the situation is much more reach and interesting come 
back to Fig. 2. Note that the perturbed timefront be- 
gins earlier than the unperturbed one and six earliest 
segments in the lower panel of Fig. 2 with identifiers 
— 111, ±112, ±113, and ±114 have no counterparts 
in the upper panel. The point is that rays with these 
identifiers in the unperturbed waveguide have launch 
angles exceeding the maximum launch angle in our 



3 Action-angle variables 

3.1 Range-independent waveguide 

First, we define the action-angle variables in a range- 
independent waveguide with c = cq(z) and, corre- 
spondingly, n — c r /ca(z) — no(z). In such a waveg- 
uide the Hamiltonian remains constant along the ray 
trajectory: 

Ho = -^n%{z)- V \ (9) 

This is the SnelFs law fl4j (in geometrical optics 
it is often presented in the form rio(z) cosx = const) 
analogous to the energy conservation law in classical 
mechanics. Equation (^ establishes a simple relation 
between the momentum p and coordinate z 

p = ±y/n 2 (z)-Hl (10) 
The action variable / is defined as the integral |9[ 
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and 



1 

2^ 



pdz = 



dz 



\l n K z ) 



H 2 



(11) 



running over the cycle of ray trajectory. Here z m - ln 
and z max are the depths of upper and lower ray turn- 
ing points, respectively. Equation ( |l"l"l) determines 
"energy", H Q , as a function of the action variable, /. 
Note the relation 



dH 2tt 



dl 



— = w, 
D 



where 



D = -2H n 



dz 



H 2 



(12) 



(13) 



is the cycle length of the ray path, and u> is the an- 
gular frequency of spatial path oscillations. Equation 
([[f) follows from (|), ©, and 

The canonical transformation from the position- 
momentum, (p, z), to the action-angle, {1,0), vari- 
ables 

I = I(p,z), 6 = 6{p,z) (14) 

and the inverse transformation 

z = z(I,6), p = p(I,6) (15) 

are determined by the equation || 

dS = pdz - H a dr = dG- 6dl - H dr, (16) 

where G = G(I, z) is the generating function. An 
explicit expression for G(I,z) is well-known J|, |ic|| . 
We represent it in the form 



G(I,z) 



J z z _ dz^n 2 {z )-H 2 {I), p > 
2nl- J z Z itiin dzy/n 2 (z)-H 2 (I), p<0 

(17) 



where z m - ln and z max are considered as functions of I. 
Then, equations 



dG 

dz 



±<Jn 2 (z)-H 2 (I), (18) 



dG 
dl 



2ttH (I) 
D 



J: 



2tt- 



n ^nl(z)-Hi(iy 
2ttH (I) rz dz 
D Jz, 



y/nl(z)-HZ(I)' 



p>0 
p<0 



(19) 



define the transformations (|14| ) and (|15|). 

Note, that the so defined angle variable varies 
from to 2-7T at each cycle of the trajectory. We 
assume that the cycle begins at the minimum of the 
trajectory. To make the angle variable continuous, 
its value should be increased by 27r at the beginning 
of each new cycle. It should be emphasized that both 
functions in Eq. ( |T^ ) are periodic in 9 with period 
2tt. 

The ray equations in the new variables take the 
trivial form 



dJ 
dr 



dH Q 

de 



0. 



de 

dr 



dH 
dl 



0,(7) (20) 



with the solution 



/ = Is 



w(I s ) r, 



(21) 



where I s and 9 S are starting values of the action and 
angle variables, respectively, at r = 0. 

Let us emphasize an almost trivial point which, 
nonetheless, is crucial for our subsequent analysis. 
Although canonical transformations ( |l^ ) and ([l5]) are 
determined by the function no(z), formally, they can 
be applied in a waveguide with a different refractive 
index profile. Moreover, these transformations can 
be used in a range-dependent waveguide, as well. 

3.2 Range-dependent waveguide 

Turn our attention to a more realistic model of the 
sound speed field given by Eq. (||). In this range- 
dependent environment with n(r, z) = c r /c(r,z) we 
define the action-angle variables using the same re- 
lations as in the unperturbed waveguide with the re- 
fractive index no(z) = c r /co(z). 
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Rewrite the Hamiltonian H 
form 



H = Hn 



V, 



= — \J n 2 — p 2 in the In what follows we shall consider 9 as a continuous 
variable defined in accordance with the remark made 
after Eq. @. Then 

0(0) = S , 9{r)=2nN + 9 e -9 s , (28) 



(22) 



and 



with Ho(p, z) defined by Eq. 

V(p,z,r) = -y/n 2 (z,r) - p 2 + ' n 2 {z) - p 2 . (23) 

If \Sc\/c r « 1, and \p\ « 1 (in underwater acous- 
tics both conditions are typically met) 



V = -Sn 7 



(24) 



where Sn = n — no — Sc/c r . Since our generat- 
ing function G does not depend explicitly on r, the 
Hamiltonian in the new variables (I, 9) is obtained by 
simple substitution of functions fll5| ) into the above 
equations. This yields 



H(I,6,r)=H (I) + V(I,9,r). 



(25) 



An explicit expression for the action function 
(eikonal) is obtained by integrating Eq. ( |l6|) with 
Hq replaced by H . We rewrite this equation as 



dS = d(G - 91) - Id9 - Hdr 



and note that even though G and 9 are discontinu- 
ous at minima of the ray trajectory, the term G — 91 
vanishes at these points and, so, this term varies con- 
tinuously along the ray. The eikonal for a ray con- 
necting points (0,z s ) and (r, z e ) can be represented 
in the form 



S = J (IdO-Hdr) 



-G(z s ,I s ) + G(z e ,I e ) + 9 S I S - 9 e I e 



(27) 



where 9 S and 9 e are angle variables at ranges and 
r, respectively. They are defined by Eq. ( p~g| ) with 
z, I, and p being ray parameters at the beginning 
and at the end of the path. Note, that the subscript 
V marks starting ray parameters while 'e' marks pa- 
rameters at the end of the trajectory. This notation 
will be used throughout the paper. 



where TV is the number of minima of the ray path, 
9 S and 9 e are the same quantities as in Eq. ( p^ ) 
(0 < 9 s ,9 e < 2tt), 61(0) and 9{r) are values of the 
continuous angle variable at the beginning and at the 
end of the ray path, respectively. 

The Hamilton equations now take the form 



dr 



and 



where 



d9 
dr 



V, 



(29) 



(30) 



89 ' 1 



dV_ 
dl' 



Two comments should be made to this definition 
of the action-angle variables. 

(i) The action variable introduced in this way (we 
]) does not conserve along 



(26) have followed Refs. (10 



the ray path even in a waveguide with very smooth 
range-dependence, i.e. our actions are not adiabatic 
invariants. Another definition of these variables (9) 
where the action does have a property of adiabatic 
invariance is shortly described in Sec. |[ 

(ii) Splitting of the Hamiltonian into a sum of the 
unperturbed constituent, Hq, and the perturbation, 
V, have been made in anticipation of our later use of a 
perturbation expansion based on smallness of 8c and, 
hence, V . However, for now we have not assumed the 
perturbation to be small and all equations derived so 
far are exact. 



4 Explicit expression for differ- 
ence in ray travel times 

In this section we compare two rays one of which 
propagates in an unperturbed range-independent 
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waveguide with Sc = 0, while another one propa- 
gates in a range-dependent waveguide with nonzero 
Sc. Both rays start at r = and our task is to de- 
rive an analytical relation for the difference in their 
travel times at a given range r > 0. Starting parame- 
ters of the rays are, generally, different but we assume 
that their action variables remain more or less close 
at any intermediate range. This assumption will be 
quantified later on. 

To distinguish between similar parameters of the 
two rays under consideration, parameters of a ray 
in the unperturbed waveguide will be marked with 
the overbar. For example, starting and final param- 
eters of its trajectory will be denoted by (p s ,z s ) and 
(p e ,z e ), respectively, while for another ray we shall 
write (p B ,z s ) and (p e ,z e ). 

The symbol A will be used to denote the difference 
between any characteristic of one ray and its coun- 
terpart for another ray, e.g. AI = 1 — 1, AS = S — S, 
Az 8 = z s — z s , and so on. 



4.1 Difference in eikonals 

Presenting Eq. ( |30| ) in the form 

d6 = (lu + Vi) dr 



(31) 



and substituting it into Eq. ( |27| ) we transform the 
expression for the eikonal to 



where 



S = So + So + Svj 

Sg = -G(z s ,I s ) 
+I s e s + G(z e ,I e )-I e 9, 

So = J {Iu - H ) dr, 



(32) 



(33) 



(34) 



S v = J (IVi - V) dr. (35) 
Equations @ and Q yield 

2ttN + 6 e -6 s = J (lu + Vi) dr. (36) 



For a ray in the unperturbed waveguide (V — 0) the 
action / does not depend on range and Eqs. (^2|)-(|36|) 
translate to 



S — Sq + S , 



(37) 



S G = -G(z s ,I) + G(z e ,I) + (6 S - 9 e )I, (38) 



S = {led- H ) dr, 



(39) 



2-kN + e -0 s = w dr. 



(40) 



Turn our attention to the difference in eikonals 
AS = S — S 



Sg — Sg + So — So + 5V- 



(41) 



Subtracting Eq. (|3§) from Eq. (|33j), exploiting Eq. 
(|Tg|), and retaining only the first order terms in AI, 
we get 

S G -S G = -G(z s , I) + G{z s J) + G(z e ,I) - G(z e , I) 



+I(A9 S -A9 e ). (42) 
Represent the difference So — So as 

So - S Q = J [F(I) - F(I)] dr, (43) 

where 

F(I) = Jw(J) - H . (44) 
Using Eq. (O), it follows that 



F(I)-F(I)=Iu ) ' + J2 -f Ar, (45) 



where 



J = W M = ^B. (46) 

dl dl y J 
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Subtracting Eq. Q from Eq. Q we get 
27rAiV + A6 e - A8 S = J [a (I) - w(J) + Vj\ dr 



W/=l , 



(47) 



Multiplying this equation by / and combining it with 
Eqs. @, @, and @ we, finally, obtain 

AS = 2nANI 



-G(z s ,I) + G(z s J) + G(z e J)-G(z e J) 




!/=2 



+ /" (A/Vr -7) dr. 



(48) 



4.2 Constituents of travel time varia- 
tion 

The difference in travel times of our two rays, At = 
AS/c r , can be represented as a sum of 4 constituents: 

At = Ate + Atjv + At i + Aty, (49) 

where 

c r At G = -G(z s , I) + G(z s ,I) + G(z e , I) - G{z e , I) 

(50) 



c r At N = 2nANI, 



(51) 



c r At i 



-u'AI 2 + -u" A/ 3 
2 3 



--w'"A/ 4 



dr. 



(52) 



c r At\ 



{AIVi - V) dr. 



(53) 



In the next section we obtain approximate solu- 
tions to stochastic ray equations in terms of action 
and angle variables. These solution provide simple 
estimations of magnitudes of individual terms present 
in Eqs. ( |52| ) and (|53|). Then the expression for At 
can be simplified by neglecting small terms. This is 
done in Sec. [| 



5 Stochastic ray dynamics in 
terms of action-angle vari- 
ables 

In this section proceeding from the Hamilton equa- 
tions ( p9[ ) and ( |3o"l ) we develop a simple approximate 
statistical description of ray trajectory fluctuations. 
This result is applied to estimate magnitudes of terms 
present in Eqs. (|5^) and ([53|) and find small parame- 
ters in the problem that may be used to simplify the 
expression for At. 

5.1 Fokker-Planck equation for the 
distribution of action 

Consider the perturbation V(I,6,r) as a random 
function with given statistical characteristics. Then 
the Hamilton equations ( |29| ) and (|o|) constitute a 
system of stochastic (Langevin) equations. We as- 
sume that the medium inhomogeneities are so weak 
that the horizontal scale of their fluctuations, l m , is 
much less than the scale of range variations of the 
action, li, 



l m "C //• 



(54) 



This assumption allows one to derive the Fokker- 
Planck equation for the probability density function 
(PDF) of the action /. Let us derive this equation 
using the method described in Ref. Jj8| . Consider an 
arbitrary smooth function (f)(1) and derive an equa- 
tion for < 4> >! where the symbol < ... > means 
statistical averaging over random medium inhomo- 
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geneities. Exploiting Eq 

d 



dr 



<(/)>- 



we get 

<<f> I (I)V e (I 1 9,r)> . (55) 



Here we have exploited the relation 



Here and in what follows the subscripts I and 9 at V 
denote partial derivatives of function V(I, 8, r) with 
respect to the first and second arguments. Transform 
the right hand side of Eq. ( |55| ) exploiting the condi- 
tion ( |54| ) . The latter ensures that there exists a range 
r < r such that 



l m < r - r < h 



(56) 



Denote the values of the angle and action variable 
at the range r$ by 9q and Iq, respectively. These 
variables at the range r may be represented in the 
form 



I = I (o) +J (1) = 0(0) +0(1) 



(57) 



where 7 (0) and 0(°) are solution of Eqs. @ and @ 
with initial conditions I^°'(ro) = Iq and 0™(ro) = #o- 
It is clear that 

/W=7 , (o) = O W(J ) (r - r ). 

The symbols I^ 1 - 1 and 0( 1 ) denote first order correc- 
tions due to the small perturbation V. By definition 
7«(r ) = and (1) (r o ) = 0. Then 

f ui = - / 1 V e { °\r')dr'. (58) 



The superscript (0) at Vg as well as at other partial 
derivatives of V and at V itself (see below) means 
that the arguments I and 9 of the corresponding 
function should be replaced with j(°) and 6^°\ re- 
spectively. For example, 

V^(r) = V(I ,e + u;'(I )(r-r ),r). 

Using this notation we get 

0«=a/(J o ) f I {1 \r')dr' + f dr 1 ' V} 0) '(/). 



After simple algebra the above expression may be 
presented in the form 

A 



0(1) = JLfl(o) 
d/ 



d r 'y(°)(r'). (59) 



Wo 



= u/ 0) 



u 



(0) 



g0(Q) 
Wo ' 



(60) 



The statistical averaging denoted by the symbol 
<> can be considered as an averaging over ray pa- 
rameters at the range ro, i.e. over Iq, 9q, and over 
medium inhomogeneities V located within the range 
interval (rp, r). Let us present the term on the right 
of Eq. ( |55| ) as 

« 4>i(I)Vg{I,e,r) >e ,v>i 

where the inner brackets denote averaging over 9q 
and the inhomogeneities V, while the outer brackets 
denote averaging over Iq. Let us first consider the 
conditional average denoted by the inner brackets. 
Using the smallness of /(^(r) and 9^(r), which are 
terms 0(V), we can simplify the conditional average 

< MI) Ve(I, 9, r) >g a y = < 0/ (/ (0) + 
x^(j( o ) + lW,# + 0<V) >e ,v 

= <f>n(Io) <I il) V g (0) >g . V 

+Mh) < V e { ?I (1) + V$9<V >g y ■ (61) 

The above expression provides an approximation to 
the right hand side of Eq. (|55| ) accurate up to terms 
0{V 2 ). Using Eq. © we get 



(l (1) y>0 o ,v. (62) 



<^>e^~\l< 



Substituting this into the last average on the right of 
Eq. ( |6l| ) and using Eq. ( |59| ) we obtain the following 
expression 



< 



yC)j(i) + vi°)0(i) > fl „-< /(i) —v {0) 
v gl i -r v gg v ^6o,v— ^ 1 v g 



+V 



(o) 



dl 



dr'V^(r') >e ,v 



(63) 
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We assume that the angle variable 9q is uniformly 
distributed over the interval (0,27r). It means that 
the statistical averaging over O is defined by 



< ... >6 



2ir 

-I ...d0„ 



This yields 



< 



d2 C ir) V (0 Hr>)> 



del 



0o,V 



for n = 0, 1, ... , Eq. (p8|) implies that W satisfies 



dr 2dl dl ' 



(69) 



This is the Fokker-Planck equation for the PDF of 
the action. A more convenient expression for the dif- 
fusivity B may be obtained by substituting Eq. dHF 
into Eq. ©. This yields 



B = 2 I dr' < V e ( °\r)V e w (r') >V 



(70) 



dV (Q) {r) dV W(r>) 
= - < — — — >e ,v (64) 



de ae 



Using relations 



vP -*;*>», vff-^vm (65) 

(note that the partial derivative with respect to I® 
has not this property: according Eq. ( |60| ) Vj ^ ^ 
^V (0) ) we obtain 

s i/(°)r(i) , i/Cl/ili) ^ _ 9 ^ r(i)T/(°) ^ 



dl 2dr V 



>B ,V 



(66) 



We assume that our perturbation is statistically 
uniform along the range r. Then the right hand sides 
of Eqs. ( |62| ) and ( |66| ) does not depend on range. 
Introducing the notation 



(67) 



combining Eqs. (|55|), (pl|), (|62|), and ( p6| ) and replac- 
ing Jo with I we arrive at 

d<<P> 1 ^ R rf 2 ^ 1 dB d(j> 
^—=2 <B dP > -2 < ^TdI > - (68) 

Consider the probability density function (PDF) of 
action W(I,r). Since <j) is an arbitrary function of I 
and 

d™6 f B n 



Since < Vg (r)Vg (r') > becomes negligible if r — 
r 1 ^> l m and the above integral does not depend on 
the lower limit r^ provided r — ro ^> l m , we can for- 
mally replace rg in Eq. ([701) with — oo. 



Let us rewrite Eq. ( ]70| ) in the form convenient for 
numerical evaluation of B. Note that for the unper- 
turbed trajectory with the action I 



dz 1 dz 
86 ~ u(j)~dr 



1 



P 



P 



>{I) ^ 



(71) 



Exploit the approximation ( p4| ) and denote by Z(I, r) 
and P(I, r) the coordinate and momentum of an un- 
perturbed ray trajectory with the action /, respec- 
tively. Then Eq. (|70| ) can be rewritten as 



B(I) 



dpK e e(I, p), 



(72) 



where 



K ee (I,p) 



Rco 2 (I) 



x / drSn z {r)Sn z (r - p))P{I,r)P{I,r - p), (73) 



with 



Sn z (r) 



ddn(r, z) 



'dz 



(74) 



z=Z(I,r) 



It is assumed that the integral should be numerically 
calculated for a typical realization of the inhomogene- 
ity Sn over a range R much greater than both the cor- 
relation length of Sn(r, z) and the cycle length of the 
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unperturbed ray. The result should not depend (at 
least, should not depend significantly) on the starting 
angle variable of the unperturbed ray and on a par- 
ticular realization of the inhomogeneity used in the 
numerical calculation. 

5.2 Statistical characteristics of per- 
turbation 

Evaluation of < V > . Consider the mean value V 
at a fixed range r. Note that even if the approxima- 
tion ( |24j ) is valid, i.e. V = —5n, and < Sn(r, z) >= 
at any fixed point (r, z) (we assume that this con- 
dition is met), generally, < V(I,0,r) >= — < 
Sn(r,z(I,9) >^ 0. The point is that I and 9 at 
the range r depend on all inhomogeneities at ranges 
r' < r, including those located within the interval 
r — l m < r' < r. It means that I and 6 are, generally, 
correlated with Sn at the range r and when averaging 
Sn(r, z(1, 9)) we cannot consider the argument z(I, 9) 
as a constant. 

In order to find the mean perturbation < 
V(1, 9, r) > we again express / and 9 at the range 
r through their values at the range r$ satisfying the 
condition (|5^) and will consider conditional averages 
for a fixed value of Iq ■ Applying the expansion ( |57| ) 
we find that up to terms 0(V 2 ) 

< V(I,9,r) > V , 0O =< V(lW,6l°\r) > v ,e 



By analogy with (p4h we have 



< 



= ^>„„, (77) 

Using this relation combined with Eq. (|65|), Eq. ( |76j ) 
may be transformed to 



<V(I Q ,9,r) > 



V,8 



d 



dloJ- 



= ir dr'<V e {0 \r)V^(r')> eo ,v 



(78) 



A further averaging over Iq provides < V >. 

In the same manner as it has been done for the 



diffusivity (see the end of Sec. 5.1) we present the 
expression for < V > in the form 



<V>= J dIW{I,r)—A 1 {I), 



where 



MI) = \I dpK g (I,p), 
^ J-00 

than for steep rays 



(79) 



(80) 



< 



>v,t 



(75) 



K e (I,p) 



Rlo(I) 



Since /W and are statistically independent of in- 
homogeneities at the range r, the first term on the 
right of Eq. ( |75] ) vanishes. The second term may be 
transformed in the same manner as we have trans- 
formed the terms on the right of Eq. (61). Using 
Eqs. (H)-@ we get an analog to Eq. (p|) 

<V(I ,9,r) yy^Kl^l-vio) 



+V {0) — 
+Ve dl 



dr'V^{r') > 9oy 



(76) 



< / dr 5n z (r)Sn(r — p)P(I 1 r). (81) 
'0 

Evaluation of < Vi >, In a similar way it can be 
shown that 



where 



:Vi>= J dIW(I,r)-^A 2 (I), 

MI) = 11 dpK 6I (I,p), 
L J-00 



(82) 



(83) 
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x / dr 5n z (r)5n z (r - p))P(I,r - p)Z I (I,r), 
Jo 



(84) 



Zj(I,r) = —Z(I,r). 



Evaluation of < (fVdr) >. It is clear that the 
mean integral 



< / V(I, 6, r) dr 1 >= / < V(I,6,r) > dr' (85) 
Jo Jo 



may be estimated using Eqs. (|79|)-(|73|). By analogy 
with the expressions for the diffusivity and < V > 
derived above we easily obtain 



< U V(I,9,r)dr') >=r I dpK(I,p), (86) 
where 



K(I,p) = -j dr 5n(r)6n(r - p)). (87) 



Evaluation of < (J Vjdr) >. Similarly, 

2 



0.05 0.1 0.15 0.05 0.1 0.15 




0.05 0.1 0.15 0.05 0.1 0.15 




Figure 6: Statistical characteristics numerically eval- 
. „, N .. , v uated over the unperturbed ray trajectory versus the 

< ( / Vi(I, 9, r) dr' ) >= r / dpKu(I, p), action /. The integrals presented in panels c, e and f 
\"' J x are evaluated over the 3000 km range. 



(88) 



where 



Ku(I,p) 



R 



x [ dr5n z (r)5n z (r ~ p))Z I (r)Z I (r - p). (89) 
Jo 
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5.3 Approximation of the action vari- 
able by a Wiener process 

Figures 6 a-f present some statistical characteris- 
tics numerically evaluated using equations derived in 
Sees. 5.1 and 5.2. All the characteristics have been 



calculated by integrating over unperturbed trajecto- 
ries placed into realizations of the perturbed waveg- 
uide. Consistent with our expectation results only 
slightly depend on particular realizations of the in- 
homogeneity used in calculations. The characteris- 
tics are considered as functions of the action variable 
determining the shape of the unperturbed trajectory. 

In Fig. 6a it is seen that the dependence of the 
diffusivity B on I is comparatively weak and we can 
approximate it by a constant 



B = 1.5 1(T 7 km. 



(90) 



Due to the smallness of perturbation the term on the 
right of Eq. ( p9| ) varies on scales small to that of the 
range variations of / and in this sense the term Vg 
can be considered as delta-correlated. Therefore we 
can present the action as 

J(r) =I s + x(r), (91) 

where I s = I{0) and x(r) is governed by equation 

dx 



dr 



£(r), 



(92) 



where £ is a white noise with 

<£>=0, < £(r)£(r') >= BS(r - r'). (93) 



Here B is a constant defined by Eq. (]90j). Equations 
(^l|)-(|93]) idealize a random increment of the action, 
x(r), as a Wiener process [|l6| with 



< x >— 0, < x(ri)x(r 2 ) >= B^(I) min(ri,r 2 ). 



(94) 



The variance of a; is a linear function of range 



Turn our attention to the angle variable. Repre- 
senting 9 as 

9(r) =9 s +Lj{I s )r + y{r), (96) 
substituting this and Eq. ( ]9l|) into Eq. ( |3C| ) we get 

dy 



dr 



u(I s + x) - u(I a ) + Vi. 



(97) 



Simplify this equation using the following two small 
parameters in the problem. 

(i) Parameter fi. It is defined as 



/' : 



1L 

I* 



'Br. 



(98) 



with I* — l^'/^lrms being a characteristic scale of 
the smooth function oj(T). The smallness of \i allows 
one to replace u{I s + x) — oj(I s ) with ui'(I s )x. It is 
important to emphasize that, typically, I s <C /* and 
the condition 



At < 1 



(99) 



does not imply that the magnitude of action fluctu- 
ations is small compared to I s . In our model typical 
values of /x at a 3000 km range are 0.1 for flat rays 
and 0.01 for steep ones. 

(ii) Small parameter /2. This parameter esti- 
mates relative contributions to y from two terms on 
the right of Eq. @: 

<u;v ld r>) 2 >^ 



< (^x(r')dr') >V2 



<(/ r W) 2 > 
|u/| (Br 3 /3) 1/2 



1/2 



(100) 



Values of < (J Q r Vidr') > x / 2 at a 3000 km range are 
shown in Fig. 6e. 

Note that according to ( pij ) Vi = —Sn z zi, where 
Sn z = dSn/dz, zj = dz/dl. An order-of-magnitude 
estimate of zj can be obtained in the following way. 
At 8 = 7r and 9 = 2n the unperturbed ray has turning 
points. For these values of 9 we have n(z) — —H, and 



E< X 2 >- 



Br. 



(95) 



n z zi 



(101) 
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where n z = dn/dz. Denote order-of-magnitude esti- 
mates of zi, n z ,Sn and Sn z by Zj, n* z , Sn* and Sn* z , 
respectively. In our model 



0.015 1/km, 5n* = 5 ■ 10~ 7 , 



Sn* z = 0.08 1/km. 



(102) 



Since a typical value of oj is about 0.15 1/km, Eqs. 
(|ToT|) and (p02[) give z} = 10. A rough estimate of 



> 



1/2 



< (/J* Vjdr') > 1 / 2 is given by < (J Q r 5n z dr') 
times zj. Using data presented in Figs. 6e and 6f it 
follows that both this rough estimate and the direct 
evaluation of < (J Q r Vjdr') >^ 2 give approximately 
the same result. Taking the value of B from Eq. ( P0|) 
and values of a/ from Fig. 15 presented below we 
find that at a 3000 km range the parameter jl varies 
within the interval of 0.1 for steep rays to 0.01 for 
flat ones. 

Due to smallness of fi and p, we can simplify Eq. 
( p7|) by neglecting the last term on the right and 
retaining only the first order in the expansion of 
u(I s + x) — uj(I s ). This yields 



dy ,, T , 



(103) 



i.e. the random component of the angle variable is 
an integral of a Wiener process 



y(r) = uj'(Is) / x(ri)dri. 
Jo 



Its variance 



a 2 e =<y 2 >=(lu') 2 Bj. 



(104) 



(105) 



Thus, we have replaced exact ray equations ( p9| ) 
and ( |30|) by remarkably sim ple approximate stochas- 
tic equations ( |92|) and ( |103| ). 

To demonstrate accuracy of Eqs. ( |95| ) and ( |105| ) 
consider a numerical example. Figure 7 graphs stan- 
dard deviations of action and angle variables, i.e. 07 
and ag , respectively, for a ray path starting at a depth 
of 0.78 km. Solid curves present results of averaging 
over a fan of 100 rays with launch angles from a nar- 
row interval centered at \ c — 7.8° (corresponding 




3000 



3000 



Figure 7: Standard deviations of the action (upper 
panel) and angle (lower panel) variables (07 and ag, 
respectively) as functions of range computed for a fan 
of 100 rays escaping a point source set at a depth of 
0.78 km. Launch angles of the rays span a narrow 
angular interval centered at 7.8°. The solid curves 
present results of averaging over the fan rays. The 
estimation 07 given by Eq. ( |95| ) is used to produce 
the dashed line in the upper panel. The dashed line 
in the lower panel is the estimation for ag given by 



Eq. (105) 



action I c — 0.06 km). Dashed curves show depen 
dencies given by Eqs. ( 



and (105). It is seen that 
the simple statistical model considered in this section 
provides reasonable predictions for standard devia- 
tions of x and y. Similar results have been obtained 
for rays starting at different launch angles. 

Relations ( |95| ) and ( |105| ) describe statistics of a 
cluster of rays with launch angles close to some fixed 
value. This result will be used below in Sec. |7|. 



5.4 Cluster of rays with a given iden- 
tifier 

Now let us focus on a cluster of different type. Con- 
sider rays leaving a point source and arriving at the 
given range r with the given identifier. Clusters of 
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this type with identifiers +124 and +160 form two 
groups of arrivals shown by thick points in Figs. 2 
and 3, respectively. 

To obtain an approximate analytical description of 
such clusters, compare two rays - one in the per- 
turbed waveguide and another in the unperturbed 
waveguide - with close but, generally, different start- 



ing action variables. Using Eqs. ( |96| ) and ( |104|) , we 
present angle variables of these rays at the range r as 

e = e s + uj{i s )r + J{i s ) f x (n)d ri 



for the ray in the perturbed waveguide and 

6 = 9 s +Lu{I)r 

for the ray in the unperturbed waveguide. The dif- 
ference 9 — 9 can be approximated by the relation 



[w(I. + x(r t )) - w(J)] dri 



u>'(I s ) (J s -/)r+ / x{ ri )dn 



(106) 



where we have omitted the constant 9 S — 9 S . At long 
ranges this constant becomes negligible compared to 
other terms whose rms magnitudes grow, on average, 
with r. Then the equation 



1 



x(ri)dri = -I s + I 



(107) 



can be idealized as a condition that singles out per- 
turbed rays whose identifiers at the range r are equal 
to the identifier of the unperturbed ray with the ac- 
tion /. Using this condition combined with Eqs. (|9l|), 
( |92|) and (|9^), one can evaluate statistical character- 
istics of action variables of rays belonging to the clus- 
ter. To illustrate this statement we consider a couple 
of examples. 

First, evaluate the probability density function 
(PDF) P(I S ), defined in a following way: P(I s )dI s 
is a probability that a stochastic ray with a starting 
action within the interval of I s to I s + dl s meets the 
condition (]107|), i.e. this ray belongs to the cluster 



defined by Eq. (107). In other words, we consider I s 
as a random variable and our task is to find its PDF. 

It is well-known that the PDF of a random variable 
(3 can be presented in the form 



P(b) =< 5{b -f})> 



(108) 



where the angular brackets denote the ensemble av- 
eraging operation. Note, that we denote any PDF 
by the same symbol P with an argument indicating 
a specific random vari able . 

Making use of Eq. (|l08|) yields 

P{I S ) =< 5 (l a -I+- I x( ri )drA >, (109) 



where the averaging goes over trajectories of the 
Wiener process x(r) defined by Eqs. ( |9^ ) and ( |9^ ) 
with the initial condition x(0) — 0. 

Using the Fo urier -representation of the 5-function 
we rewrite Eq. (10£) as 



P(I S 
The integral 



1 

2^ 



d'ye 



(110) 
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x(r\)dri 



(111) 



is a Gaussian random variable with zero mean, < 
g >= 0. T he average in the integrand on the right of 
Eq. (110) can be rewritten as [|17] 



< e 



119 



>= 



-<g"> 



(112) 



Substituting Eq. ( [11 2[ ) in Eq. ( |110| ) we arrive at 
a Gaussian integral over 7. Evaluating this integral 
and applying Eq. (Q) to find < g 2 > yields 



P(Is 



2-kBt 



exp 



3 (I a - 1) 2 
2 Br 



(113) 



The standard deviation of J s from its mean value I 
is given by 



a ls = ^Br/3. 



(114) 
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In the scope of our statistical model the difference 
in actions AI(p) =1—1 present in equations for 
AS and At derived in Sec. ||, is a Gaussian random 
function with a zero mean whose statistical charac- 
teristics are completely determined by the correlation 
function |0 



Q(j>i,Pa) =< AJ(pi) AI( P2 ) > 



(115) 



Note, that AI(p) = I s + x(p) — I. For short, we shall 
use the notation a 12 = A/(p 12 ). Applying again Eq. 
(108) we obtain the following expression for the joint 
PDF of / s , ai, and a 2 (it is assumed that < p\ 2 < 
r): 

P(I s ,ai,a 2 ) =< 8 (I s -I + g) 



xS(ai - I s + I - x(pi)) S(a 2 - I a + I - x(p 2 )) > • 

(116) 



Then 



Q(pi,P2) = J aia 2 P{I Sl ai,a 2 )dI s daida 2 



< (x(pi) - g){x{p 2 ) - g) > 



B I - - pi - p 2 + 1 ^ 2 + min(pi, p 2/ 



(117) 



Later on we shall use this equation to calculate sta- 
tistical moments of AI in order to estimate individual 



terms in the expansion (52). In so doing we shall ap- 
ply the known properties of statistical momenta of 
Gaussian random variables pjj], In particular, we 
shall use the relation 



< AI 2 ( Pl )AI 2 ( P2 ) >=< AI 2 { Pl ) >< AI 2 { P2 ) 



> 



+2 < AI{ Pl )AI( P2 ) > 2 



Q(pi,Pi)Q{ P2 ,p 2 ) + 2Q 2 ( Pl ,p 2 ). (118) 



5.5 Stochastic instability of ray iden- 
tifier 



In Sec. |2.3| we have already seen that in the perturbed 
waveguide ray trajectories exhibit extreme sensitivity 
to starting launch angles. An impressive demonstra- 
tion of this instability is presented in Fig. 5. The 
points representing dependence of the ray identifier 
on the launch angle are randomly scattered and this 
fact suggests that rays with very close launch angles 
are practically uncorrelated and should be described 
statistically. 

Our stochastic ray theory provides a tool for quan- 
titative description of stochastic ray instability. In 
particular, it allows one to estimate chaotic spread of 
ray identifiers shown in Fig. 5. Consider unperturbed 
and perturbed rays starting at the same launch an- 
gle, Xs, and, consequently, with the same starting 
action, I s . Denoting the numbers of turning points 
of these paths by J and J, respectively, we note that 
at long ranges both J and J are large and the same 
is true of their rms difference. Then we can approxi- 
mate A J = J - J by y/2n. This yields < A J >= 
and (Taj = (< A J 2 >) x / 2 = ag/2n with ag given by 
Eq. ( |l05| ). Assuming that for most rays with start- 
ing angles close to Xs the value of J lies within the 
interval 



J - 2r7 A j < J < J + 2<JAj 



(119) 



we find an estimate for the spread of ray identifiers. 
Figure 8 presents the number of ray turning points, 
J, against the launch angle, Xs, at 1500 km and 3000 
km. For 3000 km range we have the same plot as in 
Fig. 5, but the solid line representing the dependence 
of J on Xs is slightly smoothed. Th e dashed lines 
represent limits d efine d by Eq. (11£). This result 
confirms that Eq. dll9j ) gives a reasonable estimation 
of the spread of ray identifiers. 

Stochastic dependence of the identifier on the 
launch angle can be considered from a different view- 
point. Let us fix some ray identifier and study statis- 
tics of starting actions, I s , of perturbed rays arriving 
at the given range r with this identifier. The prob- 
ability density function of I s is given by Eq. (113). 



The mean value of I s is equal to /, action of an unper- 
turbed ray which has the given identifier at the range 
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X s (deg) 



Figure 8: Ray identifier, J, as a function of the launch 
angle, \s, at the ranges 1500 and 3000 km. The 
points depicting rays at 3000 km range are the same 
as in Figs. 2 and 3. The solid lines are smoothed 
functions J(Xs) for the unperturbed waveguide. The 
dashed lines show the limits established by Eq. (|ll 



r. According to Eq. (113) I s is a Gaussian random 
variable and it is natural to expect that most part of 
rays with the given identifier at range r have starting 
action within the interval 



I - 2<J Is < I s < I + 2(7r 



(120) 



where oj s is determined by Eq. ( |114| ). Figure 9 
demonstrates that this prediction agrees with results 
of our numerical simulation. The solid lines present 
dependencies of the starting action of the unper- 
turbed ray on the number of ray turning points at 
1500 km and 3000 km, while the dashed lines indi- 
cate the borders of intervals defined by Eq. (12C). 



Consistent with our expectation, most points depict- 
ing parameters I s of perturbed rays starting upward 
against numbers of their turning points at 1500 km 
(upper panel) and 3000 km (lower panel) lie within 
areas embraced by the dashed lines. 



Figure 9: Starting action, I s , as a function of ray 
identifier, J, for the ranges 1500 km (left) and 3000 
km (right). The solid lines are smoothed functions 
I S (J) for the unperturbed waveguide. The dashed 



lines show the limits established by Eq. (12C) 



H) 6 



Small parameters in the prob- 
lem and approximate formula 
for travel time variations 



The approximate stochastic ray theory developed in 
the preceding section allows one to deduce simple an- 
alytical estimations of terms present on the right side 
of Eq. (49). Our main concern is with terms Ati 
and Aty which are given by integrals over r and, 
in principle, may be very large at long ranges. The 
objective of the present section is to compare these 
constituents of At and simplify Eq. ( |49| ) by neglect- 
ing small terms. 

(i) Estimation of At/. Divide Atj into a sum 



Atr = At 



(2) 



At 



(3) 



A< 



(4) 



where 



Atf = ^- I A/ J ,/,-. 



Atf= £. / A/ :! dr. 



(121) 



(122) 



(123) 
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A4 4) =^- /A/ 1 ,/,, 

8c r 



(124) 



(2") (3) 

and estimate mean values of components At} , At } , 



and Atf ] . Applying Eq. ( p7| ) to find < A/ 2 > and 
making use of the relation 



< AI >= 3 < AI >" 



(125) 



valid for a Gaussian rand om variable with a zero 
mean (it follows from Eq. ( |l!8|) ), yields 



< Atf > 



12c r 



-Br z 



< Atf } >= 0, < Atf > > 



(4) 



80 c, 



(126) 



£? r . (127) 



(ii) Estimation of At v. Due to smallness of AI 
we can transform the integrand on the right of Eq. 
©as 

AI Vi(I, 0, r) - V ~ -V(I, 0, r) = 5n(r, z(I, 6)). 

(128) 



Then 



At v = / 5n(r,z(I,6))dr. (129) 



The integration in Eq. (129) goes along the unper- 
turbed ray path determined by the action I. Since 
0, generally, differs from 9, this integral differs from 
that present in Eq. ( |86| ) which is evaluated over the 
unperturbed ray trajectory. However it is clear, that 
the statistics of Atv determined by Eq. (129) is ap- 
proximately described by Eqs. ([861), and (|S7|) . Figure 
Fig. 6c shows the standard deviation of Atv/c r at 
the 3000 km range. Note that the quantity similar 
to Aty evaluated over an unperturbed ray path, is 
widely used to estimate travel time shifts due to in- 
homogeneities at sh ort en ough ranges Jj]]. 

Using Eqs. ( fl26b a nd ( fl27b we see that the first 
term in the sum ( 121 ) dominates if the parameter 



Mi 



< Atf > 



< Atf 



> 
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in 



to 



(130) 



is small compared to unity. Taking into account Eq. 
d99) we can roughly estimate /ii as 



(131) 



where /** is some characteristic scale. It turns out 
that for our environmental model /** differs consid- 
erably from the scale I* present in Eq. (p^). There- 
fore, the parameter /ii cannot be estimated as /i 2 as 
it might be expected. Nevertheless, smallness of both 
fi and Hi is caused by the same factors: the weakness 
of the perturbation combined with the smoothness of 
the function w(/). 

The smallness of \x\ allows one to replace the con- 
stituent of the difference in ray travel times Atj (see 
Eqs. (||) and (||)) with Atf\ Taking into account 
(129) we can now represent the general expression for 
the difference in travel time (49) in the form 



c r At = 2tt ATV / + c r At G + y J AI 2 dr + J Sn dr. 

(132) 

This formula is the main result of the present work. 
Note that it can be further simplified if we consider 
two eigenrays coinciding the same pair of the end- 
points in the perturbed and unperturbed waveguides. 
In this case z s = z s and z e = z e and according Eq. 



Ate = 0. 



(133) 



z a , and arrive 



If both rays leave the same point, z, 
at two different points whose depths are close, then 
using Eq. (na) we can approximately present A£g as 



At G 

If parameter 

< (At v ) 2 > 1/2 



< Atf > 



PeAz e 



12 Cr < (Aty) 2 >V 2 

Br 2 



(134) 



(135) 



is small compared to unity, then the last term on the 
right of Eq. (132), Aty, can be neglected because it 
is small compared to At/. 
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Then the fourth term on the right in Eq. (132) van- 
ishes, the action I (like I) does not depend on r, 
and the same is true of AI. We shall compare rays 
arriving at the same points, i.e. eigenrays, and in ac- 
cordance with Eq. ( po| ) neglect the term Ate- This 
yields 



c r At = 2irl AN - 



-AFr. 



(136) 



In the case AiV — N — N = 1 this equation provides a 
difference in travel times of two eigenrays with iden- 
tifiers ±J and ±(J — 2). According to Eqs. ( |2l] ) and 
( |2S| ) the action variable of these eigenrays satisfies 

- w(J) = ui' AI + 0(AI 2 ) 



Figure 10: Dimensionless parameters defined by Eqs. 



(13C) and (135) as functions of launch angle at the 
3000 km range. 



Figure 10 shows parameters /ziand [ii at the range 
3000 km as functions of the launch angle for a point 
source set at a depth of 0.78 km. It is seen that while 
fii is everywhere small, the parameter ^2 is small only 
for flat enough rays. 

In what follows we shall apply Eq. (132) to 



-(2tt- A9 s + A6 e 



(137) 



study variations of the timefront due to internal- 
wave- induced random inhomogeneities. Our primary 
concern will be with the estimation of widening and 
bias of the timefront segments in the presence of per- 
turbation. 



Analytical description 
timefront structure 



of 



7.1 Timefront in a range-independent 
waveguide 



Although Eq. (132) has been derived to compare 
travel times of perturbed and unperturbed rays, it 
can also be applied in the case when both rays prop- 
agate in the unperturbed waveguide with Sn = 0. 



Assuming that ui as a function of I is monotonous at 
the interval (/,/), at long ranges (N » 1) we have 
|A0 S | , |A0 e | < 2tt and 



to r \ r z 



Then Eq. (^36J) reduces to 

c r At = 7r(J + I). 



(138) 



(139) 



A similar result have been obtained in Refs. [|l9| 
|cj |lj] (in Refs. Q U it has been derived for an 
adiabatically range-dependent waveguide). 

An interesting a nd s omewhat surprising fact fol- 
lowing from Eq. ( ]l3S| ) is that there exists a con- 
servation law for temporal shifts between timefront 
segments. Consider a bunch of rays with launch an- 
gles within a narrow interval. Action variables of all 
these rays are close to some value which we denote by 
Jo- Beginning from a certain range r* these rays will 
form at least two segments with identifiers that differ 
by 2. When estimating the temporal shift between 
two such segments we shall compare rays arriving at 
the same depths. With this in mind, we can estimate 
the temporal shift as tq — 2irIo/c r . It should be em- 
phasized that tq does not depend on range. It means 
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that although the number of segments formed by rays 
with launch angles from a given narrow angular in- 
terval grows linearly with range, temporal shifts be- 
tween neighboring segments (to be more precise, be- 
tween segments corresponding to identifiers ± J and 
±(J — 2) with J depending on range ) remain approx- 
imately the same at any distance. A more detailed 
discussion of this issue is given in Refs. [[l9[ 

The above statement means that a simple evalu- 
ation of the action variable as a function of launch 
angle, /(x), gives a considerable quantitative infor- 
mation on temporal structure of the pulse signal valid 
at arbitrary (long enough) range. If N ^> 1 , the val- 
ues of / and I on the right side of Eq. (139) are close 
and this equation can be approximated by 



At 



Xs +Xs 



(140) 



where \ s and \s are launch angles of rays under con- 
sideration and 



-(Xs) = 2ttI(xs)/cv 



(141) 



A solid line in Fig. 11 graphs t(\s) for our model of 
range-independent waveguide. Looking at this curve 
we can predict that, for example, a difference in travel 
times of two eigenrays with launch angles close to 7° 
and numbers of cycle which differ by 1, will be close to 
r = 0.19 s at any range and at any depth, provided 
such eigenrays arrive at the observation point. 

Let us select some reference depth, z r , and define 
the temporal shift between segments with identifiers 
± J and ±( J— 2) - we denote this shift by I±j j ±(j_ 2 ) 
- as a difference in travel times of two eigenrays with 
these identifiers arriving at the depth z r . Arrivals 
of these rays in the upper panels of Fig. 2 and 3 
can be found as intersections of the corresponding 
segm ents with a horizontal line z = z r . According to 
Eq. (141) the value of T±j.±(j-2) does not depend on 
a particular value of z r (the only requirement is that 
both segments must intersect the line z = z r ). This 
result agrees with the fact that neighboring segments 
in Figs. 2 and 3 with inclinations of the same sign 
are almost parallel. 

Since T±j.±(j_2) represent a difference in travel 
times of two eigenrays, it can be estimated as 




Figure 11: The solid curve is r = 2TrI(x)/c r from 
Eq. (141). The circles and triangles are time delays 
between segments with the identifiers +J and — 
2) at a reference depth of 0.78 km computed for even 
J at 1500 km and 3000 km ranges, respectively. 
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T (x±J,±(J-2)) with X±J,±(J-2) being a half-sum of 
launch angles of the corresponding eigenrays. It 
means that points depicting values of T±j,±(j-2) 
against x±j,±(j-2) should lie on the curve t(x) at 
any range and for any reference depth, z r . In Fig. 11 
this prediction is verified for the timefront shown in 
the upper panels of Figs. 2 and 3 and for a similar 
timefront at the range of 1500 km. Circles and tri- 
angles depict 7 1 + j i+ (j_ 2 ) as functions of x+J,+(J-2) 
at 1500 km and 3000 km ranges, respectively, for 
z r = 0.78 km. Travel times shifts for only even J 
(from 60 to 96 at 1500 km, and from 118 to 196 at 
3000 km) are shown. It is clearly seen that all the 
circles and triangles are, indeed, located close to the 
solid curve representing r(x). 

7.2 Timefront in the presence of per- 
turbation 

7.2.1 Widening and bias of timefront seg- 
ments 



In the presence of weak range-dependent inhomo- 
geneities the structure of timefront becomes more 
complicated: instead of infinitely thin segments of 
smooth curves, we have some areas filled with ran- 
domly scattered points. Although we observe the 
scattered points only because our fan is far too sparse 
to resolve what should be unbroken curves, the ap- 
pearance of such regions indicates the presence of 
chaotic rays. 

As it has been pointed out in Sec. 



2.3 the early 



portion of the timefront formed by steep rays still 
"remembers" its structure in the unperturbed waveg- 
uide. The points depicting arrivals of rays with the 
given identifier are scattered in the vicinity of the cor- 
responding unperturbed segment. A group of arrivals 
formed by rays with the same identifier produces a 
fuzzy versions of an unperturbed segment. We shall 
call such groups of points in the time-depth plane, 
the fuzzy segments. Thus, every fuzzy segment, like 
every segment of the unperturbed timefront, is asso- 
ciated with some identifier. Two examples of fuzzy 
segments are shown in Figs. 2 and 3. We mean 
two groups of thick points indicating arrivals of rays 
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Figure 12: Arrivals with identifiers +124 (left panel) 
and +160 (right panel) at the range 3000 km pre- 
sented in the time-depth plane. The points and solid 
lines depict arrivals with and without internal waves 
present, respectively. The magnitude of time delay, 
|At|, between arrivals of the perturbed ray and an 
unperturbed ray with the same identifier is shown 
for the earliest arrival with identifier +160. 



with the identifiers +124 and +160. Closer views of 
these groups of points are shown in Fig. 12. The 
left (right) panel presents arrivals with the identifier 
+124 (+160). Thick solid lines in both panels show 
the unperturbed segments with the identifiers +124 
(left panel) and +160 (right panel). 

In order to derive quantitative characteristics of the 
fuzzy segment describing its spread and bias, we in- 
troduce the quantity At defined as follows. Consider 
a particular ray (perturbed) contributing to the fuzzy 
segment and denote its travel time by t p . A travel 
time of an unperturbed ray with the same identifier 
and the same arrival depth we denote by t u . Then 

Ar = t p ~t u . (142) 

In words, At represents the distance along the i-axis 
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Figure 13: Time delays, Ar, defined by Eq. ( |142| ), 
for arrivals of rays with identifiers +124 (left panel) 
and +160 (right panel) at the range 3000 km. The 
solid lines represent results of the direct ray tracing, 
while the dashed lines are predictions made using Eq. 
(143). For rays with identifier +160 only the term 
At/ in Eq. (|143|) has been taken into account. 



between the given point of the fuzzy segment and the 
unperturbed segment with the same identifier. In the 
right panel of Fig. 12 the magnitude of At is shown 
for the earliest arrival. Note, that Ar is defined only 
for those rays forming the fuzzy segment whose ar- 
rival depths lie within a depth interval covered by the 
unperturbed segment. An approximate analytic al ex - 
pression for Ar is readily obtained from Eq. (132). 
Since we consider eigenrays with identica l ide ntifiers, 
the first two terms on the right of Eq. (132) vanish 
and we arrive at 



where 



and 



Ar = Ar/ + At v , 



Ar 7 = ^— / A/ 2 dr, 



(143) 
(144) 



Ar y = — / 5n(r,z(r))dr. (145) 

Cr J 



First of all, check an accuracy of Eq. ( |143| ) at 3000 
km range. This is done in Fig. 13 where we present 
travel time shifts Ar for arrivals shown in Fig. 12. In 
the perturbed waveguide there are 160 fan rays with 
the identifier +124 whose depths at 3000 km belong 
to the depth interval covered by the unperturbed fan 
rays, i.e. between upper and lower points of a solid 
curve in the left panel of Fig. 12. A solid line in the 
left panel of Fig. 13 connects exact values of Ar ob- 
tained by ray tracing, while a das hed line represents 
predictions provided by Eq. (143). Launch angles of 
perturbed rays with the identifier +124 belong to the 
interval (7°, 9.5°). In Fig. 10 we see that the param- 
eter ^2 for such launch angles cannot be considered 
as small comp ared to unity, which means that both 
terms in Eq. ( |l43j ) should be retained. In the right 
panel of Fig. 13 a similar plot is shown for 139 rays 
with the identifier +160. The prediction depicted by 
a dashed curve has been made by retaining only the 
term Ar/, because launch angles of these rays are 
less than 7°. The parameter /j.2 for such launch an- 
gles is small (see Fig. 10) and the term Ary can be 
neglected. Figure 13 demonstrates that Eq. (143) 
provides a reasonable estimation for Ar. 

Figure 12 exhibits a new phenomenon. The per- 
turbation causes not only diffusion of the timefront 
segment but it also leads to some regular bias: rays 
with the given identifier in the perturbed waveguide 
arrive, on average, earlier than unperturbed rays with 
the same identifier. A similar bias is observed for 
every fuzzy segment. A qualitative explanation to 
this effect follows i mme diately from the fact that the 
term Ar/ in Eq. ( |l43|) usually dominates. This is 
true even of steep rays although in this case the term 
Ary should be retained for obtaining an accurate pre- 
diction of Ar. The sign of Ar/ is determined by the 
sign of the derivative w' . The latter is negative for 
rays propagating without reflection off the surface 
and bottom, because in typical deep ocean waveg- 
uides the cycle length of the refracted ray grows with 
the launch angle. The spatial frequency uo and its 
derivative with respect to / are shown in Fig. 14 
for the model of unperturbed waveguide on which we 
rely in this paper. 

Approximating Ar by Ar/ (i.e., practically, by 
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Figure 14: Left panel: Spatial frequency of ray tra- 
jectory oscillation, oj, in the unperturbed waveguide 
as a function of the action variable, I. Right panel: 
The derivative of the spatial frequency uj with respect 
to /. 

At) '), one can estimate a mean bias of the f uzzy 
segment, < At >, by making use of Eq. ( |l26| ). In 
Fig. 15 we compare this prediction (solid line) with 
"real" mean values of At (circles) evaluated for fuzzy 
segments formed by our fan rays with identifiers + J 
at ranges of 1500 km and 3000 km. 

The rms time spread of the fuzzy segment can be 
estimated as 



y/< (At- < At >) 2 > 



V<A 



T > 



< At > 2 



Appr oxim ating again At by At/ and making use of 
Eq. ( |ill) we get 

' Ar 2 >= (^-] j / <!,,,<!<>, 



o Jo 



x [Q(pi,pi)Q(p 2 ,p 2 ) + 2Q 2 (pi,p 2 )] = 
This yields 

Iw'l Br 2 . . 

ctat = 1 _ ,- = 0.9 | < At > 



2c, 
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Figure 15: Mean bias of fuzzy segments due to in- 
homogeneities as a function of the ray identifier at 
ranges of 1500 km (upper panel) and 3000 km (lower 
panel). Each circle shows the bias averaged over a 
fuzzy segment, i.e. over a group of rays with the 
given identifier + J. Th e solid lines are predictions 
obtained using Eq. (126). 



i.e. the rms time spread is prac tical ly equal to its 
mean bias. In other words, Eq. ( |126| ) estimates not 
only the bias of the fuzzy segment but its time spread 
as well. In Fig. 16 we compare this prediction to 
"real" time spreads obtained in numerical simulation. 
Results presented in Figs. 15 and 16 demonstrate 
that even though predictions made with Eq. (126) 



have some systematic error for very flat rays (corre- 
sponding to large J in the lower panels of Figs. 15 
and 16), these predictions are in reasonable agree- 
ment with numerical simulation and give the correct 
order of magnitude for both time spread and bias. 
In particular, we see that at the 3000 km range, in 
accordance with Eq. (126), the bias and time spread 
become 4 times larger compared to their values at 
1500 km. 

Looking at Figs. 15 and 16 we see that the bias and 
time spread are especially small for segments formed 
by steep rays (small J). On the right hand side of 
Eq. ( |l26j ) there is a factor uj' , depending on I and, 
hence, on the launch angle. Looking at Fig. 6 and 
at the right panel in Fig. 14 we conclude that the 
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dependencies of < At > and ctat on the launch angle 
are mainly determined by the factor uj'. 
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Figure 16: Width of a fuzzy segment as a function of 
the ray identifier at ranges of 1500 km (upper panel) 
and 3000 km (lower panel). Each circle shows the 
standard deviation of the time delay t obtained by 
averaging over the fuzzy segment, i.e. over a group 
of rays with the given identifier + J. Th e solid line s 
are predictions obtained using Eqs. ( |l46| ) and (126). 



7.2.2 Resolution of fuzzy segments in the 
perturbed timefront 

In the upper panels of Figs. 2 and 3 it is clearly 
seen that unperturbed segments come in groups of 
four and each group consists of segments with the 
identifiers — (J — 1), ±J, and +(J— 1). The function 
t(x) defined by Eq. (141) predicts the time delay 



between two consecutive groups of four formed by 
rays with launch angles close to x- The difference 
in travel times of two neighboring segments can be 
roughly estimated as r/4. Estimating the width of 
the fuzzy segment as 3cr T , we introduce the parameter 



R = 



12cr T 



24tt/ 
\u>'\Br 2 



(147) 



representing the ratio of the segment width to the 
time delay between neighboring segments. If R > 1, 
the fuzzy segment is resolved in the perturbed time- 
front while the segment with R < 1 overlaps with 
neighboring segments. From the viewpoint of eigen- 
rays, the condition R > 1 (R < 1) means that 
neighboring unperturbed eigenrays in the presence of 
perturbation split into nonoverlapping (overlapping) 
clusters. 

In order to verify this prediction, consider the ray 
travel time, t, in the unperturbed waveguide as a 
function of the launch angle, Xs- For launch angles 
of the same sign the function t(x s ) is monotonic and 
the i nverse function Xs{t) is unambiguous. Equation 
( 147 ) defines the parameter R as a function of the 
launch angle Xs- Replacing Xs with Xs(t) we obtain 
the function R(t) that associates a value of parameter 
R with every segment of the unperturbed waveguide. 
If this value is greater than unity we expect that the 
corresponding fuzzy segment will be resolved. This 
statement is illustrated in Figs. 17 and 18. The func- 
tions R(t) at r = 1500 km and r = 3000 km evaluated 
for rays starting upward (for rays starting downward 
the results are practically the same) are plotted in the 
lower panels of Figs. 17 and 18. These functions are 
shown at time intervals where R(t) passes through 
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Figure 17: Upper panel: Fragment of the timefront 
at 1500 km range where transition from resolved to 
unresolved (overlapping) fuzzy segments is obse rved. 
Lower panel: Parameter R defined by Eq. (147) as a 
function of travel time. Theory predicts that resolved 
(unresolved) fuzzy segments are located to the left 
(right) of the travel time at which R passes through 
1. 



1. The corresponding portion of the perturbed time- 
fronts from the upper panels of Figs. 2 and 3 are 
shown in the upper panels of Figs. 17 and 18. It 
is clearly seen that, indeed, the condition R = 1 al- 
lows one to estimate the critical travel time which 
divides the perturbed timefront into parts consisting 
of resolved and unresolved fuzzy segments. 

The fact that early arriving fuzzy segments formed 
by steep rays are well resolved is linked to the follow- 
ing two factors. First, the function t(x) grows with 
\x\- It means that the delay between two consecu- 
tive segments of the timefront formed by steep rays 
is larger than that for consecutive segments formed 
by flat rays. Second, the rms width of the fuzzy seg- 
ment, a at, is especially small for steep rays (small J 
's), as it is seen in Fig. 16. 




Figure 18: The same as in Fig. 17 for 3000 km range. 

7.2.3 Stability of fuzzy segments and Fer- 
mat's principle 

Probably, the most surprising feature of the per- 
turbed timefront is unexpectedly small widenings of 
timefront segments. In particular, in the upper panel 
of Fig. 16 we see that at 3000 km range the maximum 
rms widening of the fuzzy segment does not exceed 
0.12 s. On the other hand, Fig. 4 demonstrates that 
a typical time spread for a cluster of rays with close 
launch angles is about 2 s, i.e. much larger. Loosely, 
we can state that the ray travel time dependence on 
the ray identifier is less chaotic and much more pre- 
dictable than its dependence on the launch angle. 

In order to inte rpret this phenomenon, we come 
back to Eq. (143) and recall that the first term on 
the right side dominates. It means that the travel 
time shift between perturbed (P) and unperturbed 
(U) rays connecting the same points and having the 
same identifier can be approximately written as 



Ar = -[So(P) - S (U)} 

C r 



(148) 



where So(P) and Sa(U) are the values of the func- 
tional So — j{pdz — Hodr) evaluated over the tra- 
jectories of perturbed and unperturbed rays, respec- 
tively. According to Fermat's (Hamilton's) principle 
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H [L2], [L4|, the unperturbed ray provides a station- 
ary path of the functional Sq, This fact explains the 
absence of the linear in AI term in Eqs. (p^), (132), 
and Eq. (|143|) . Since A/ is our small parameter, the 



absence of terms O(AI) gives some qualitative inter- 
pretation of smallness of At. In this sense, the small 
time spread of clusters of rays with the same identi- 
fier can be interpreted as a consequence of Fermat's 
principle. 

Note that the difference in travel times of rays with 
different identifiers (N ^ 0) is defined mainly by the 
term At^ given by Eq. ( |5l|) which usually is signifi- 
cantly larger than Aij. 

8 Action-angle variables in 
a waveguide with a range- 
dependent background sound- 
speed profile 

An ocean-acoustic propagation model with the 
sound speed field being a superposition of a 
range-independent background and a weak range- 
dependent perturbation responsible for emergence of 
ray chaos may be too idealized. In this section we 
shortly outline how the results obtained in the pre- 
ceding sections can be generalized to a more realistic 
model. 

First, let us shortly discuss a method of introducing 
of the action-angle variables in the range-dependent 
waveguide J|] without dividing the Hamiltonian into 
a sum of an unperturbed term and a perturbation. 
Define canonical transformations ( |l^ ) and ( |l5| ) at a 
current range r using Eqs. ( |l8| ) and (|l9j) evaluated 
for an auxiliary range-independent waveguide with 
the same cross-section that the real waveguide has at 
the range r. In this case the canonical transformation 
will be different at different ranges and Eqs. ( |l4| ) and 
( |l5| ) translate to 



6{p, z, r) 



and 



z = z(I,6,r), p = p(I,6,r). 



(149) 



(150) 



The generating function G now becomes a function 
of not o nly I an d z, but of r, as well. However, 
H = — \Jn 2 — p 2 in the new variables is a function of 
I and r, but not ||. 

The Hamilton equations in the new variables pre- 
serve their canonical form 



dl 

dr 



dH s 



d6 dH s 



89 ' dr 



dl 



with the new Hamiltonian 



H s (I,e,r)=H(I,r) + A(I,e,r) 



where 



A(I,9,r) 



dG(I,z,r) 



dr 



(151) 



(152) 



(153) 



z=z(I,6,r) 



The term A is small and can be neglected if range 
variations in the environment are adiabatic, i.e. if 
variations in the environment are small at the cycle 
of the ray trajectory. Then, dl/dr — and I re- 
mains constant along the ray trajectory, i.e. the ac- 
tion variable defined in this way does have a property 
of adiabatic invariance. 

However, we suppose that this common approach is 
not convenient for description of the chaotic ray mo- 
tion induced by random internal waves. The point 
is that if A is not negligible, then the connection be- 
tween H, A and 5c becomes n on-tr ivial, and it is dif- 
ficult to divide Hamiltonian fll52j ) into a sum of a 
smooth unperturbed term and a small perturbation. 
But such a decomposition of the Hamiltonian is nec- 
essary for application of our perturbation theory. 

A more appropriate approach can be developed if 
the sound speed field is a sum of a smooth range- 
dependent sound speed, co(r, z), and a weak pertur- 
bation, Sc(r, z), i.e. 



c(r, z) = co(r, z) + Sc(r, z). 



(154) 



Instead of the range-independent unperturbed 
waveguide considered in the preceding sections, now 
we have an adiabatic one. Then, it is convenient to 
introduce the action-angle variables at every range r 
using an auxiliary range-independent waveguide with 
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the cross-section coinciding with that of the unper- 
turbed waveguide. This yields the new Hamiltonian 
in the form 



H s =H (I,r)+V(I,9,r), 



(155) 



where V(1, 6, r) is the perturbation defined in Eq. 
( p3| ) with the unperturbed refractive index no now 
depending not only on z but on r, as well. Then 
the Hamilton equations have the same form as Eqs. 
( p9|) and (^) in Sec. |3.2| , although the angular fre- 
quency uj now depends on r, uj(I,r) — dHo(I, r)/dl. 
All expressions for differences in ray travel times of 
perturbed and unperturbed rays derived in Sec. |] 
remain valid for this more realistic model. 



9 Summary and conclusion 



In this paper we have derived the approximate an- 
alytical approach for description of ray travel times 
and other parameters of the ray structure in deep 
ocean environment. Our results remain valid at 
ranges up to, at least, a few thousand km. The 
approach is based on the assumptions that (i) the 
perturbation giving rise to the chaotic ray motion is 
small and (ii) even at long ranges rms variations of 
the action variable are small compared to the char- 
acteristic scale of function The dimensionless 
s mal l para mete rs in the problem are given by Eqs. 
(|l30|) and (|l35|) . 



The exact expression for the difference in travel 
times of perturbed and unperturbed rays, At, de- 
termined by Eqs. ( ^9| ) - ( |5^ ) has been significantly 
simplified by expanding it in a power series in AI 
and neglecting small terms. The smallness of fluc- 
tuations of the action variable has also been used to 
simplify the stochastic ray (Hamilton) equations. It 
has been shown that the fluctuating components of 
the action and angle variables can be idealized as a 
Wiener process, and an integral of the Wiener pro- 
cess, respectively. This result, which allows one to 
evaluate (approximately) practically any statistical 



characteristic of the ray trajectory, in the present pa- 
per has been combined with an approximate formula 
for At and applied for investigation of range varia- 
tions of ray travel times. 

Our primary concern has been with the range vari- 
ations of the timefront representing ray arrivals in the 
time-depth plane. The unperturbed timefront con- 
sists of segments of smooth curves. Each segments is 
formed by rays with the same identifier. In the pres- 
ence of perturbation segments become fuzzy: arrivals 
with the given identifier form a set of points randomly 
scattered around the unperturbed segment. The time 
spread of these points turns out to be unexpectedly 
small. It is much less than a time spread of arrivals 
with launch angles within a narrow angular interval 
corresponding to launch angles of rays forming an 
unperturbed segment. The most apparent manifes- 
tation of this phenomenon is a surprising stability of 
early portions of the timefront formed by steep rays 

33 31- 

Our approach provides a quantitative description 
of fuzzy segments. It gives estimations of their widths 
and biases. Using these estimations, it follows that 
the sign of the bias is determined by the sign of 
div/dl, the derivative of the spatial frequency of ray 
oscillation in the unperturbed waveguide with respect 
to the action variable. In typical deep ocean waveg- 
uides dio/dl < for all refracted rays, which means 
that fuzzy segments have negative bias, i.e. per- 
turbed rays, on average, arrive earlier compared to 
unperturbed rays with the same identifier. It has 
been shown that the surprising stability of fuzzy seg- 
ments with respect to the perturbation is related to 
the Fermat's (Hamilton's) principle. 

The estimations derived for the timefront segments 
can be applied to study characteristics of chaotic 
eigenrays. In Ref. (see also Ref. |^|) it has been 
discovered numerically, that in the presence of a weak 
perturbation the unperturbed eigenray splits into a 
cluster of new eigenrays with close arrival times. The 
time spread of such a cluster can be estimated as a 
time spread of the corresponding fuzzy segment at the 
corresponding depth. So, our criterion of nonover- 



lapping of neighboring segments (147) provides the 
criterion of resolution for clusters of eigenrays. 
It should be emphasized that results of the present 
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work have been obtained for the perturbation induced 
by internal waves. Their generalization to the case of 
other inhomogeneities, e.g. the mesoscale inhomo- 
geneities, requires a further investigation. 
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A Appendix. Evaluation of ac- 
tion and angle variables using 
a ray code 

Explicit analytical relations connecting the position- 
momentum and action-angle variables can be derived 
only for a few simple refraction index profiles uq{z) 
ijiol , In underwater acoustics the researcher is 

usually forced to deal with refraction index profiles 
defined by spline approximations of data measured 
at discrete depths. In this case ray trajectories are 
computed using standard ray codes which allow one 
evaluate such parameters of the ray as the coordi- 
nate, momentum, and travel time [0. Fortunately, 
the same codes can be applied for evaluation of the 
action and angle variables. It will be a simpler and 
more convenient procedure compared to the direct 
application of Eqs. ( |l8| ) and (|l9|). 

As it has been pointed out in Sees. ^2 and || 
the transformation from the position-momentum to 
the action-angle variables is always defined using an 
auxiliary range-independent waveguide. So what we 
need for establishing relations between the two types 
of variables in any waveguide, is a code for computing 
the ray path in the range-independent waveguide. 

The action variable for the given trajectory can be 



easily found by evaluating the integral 

p\r) 
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o yjnl(z(r)) -p 2 (r) 



dr 



(156) 



where the functions z{r) and p(r) are the coordi- 
nate and the momentum, respectively, as functions 
of range, and D is the cycle length. Here we have 
applied Eqs. ( |ll|) and Computing ray paths cor- 
responding to different values of the Hamiltonian (||) 
and using a spline approximation we find functions 
H {I), I(H ), and w(J) = 2n/D{I). 

Then the transformation from (p, z) to {I, 9) is per- 
formed as follows. Substituting the given p and z into 
Eq. (g) yields Hq. Then, making use of the function 
I {Ho) we get the action variable, I. In order to find 
the angle variable, 9, consider a ray trajectory with 
Ho determined by the given p and z, which begins 
at its upper turning point, i.e. at z{0) = z min (recall 
that we are using z-axis directed downward) . Accord- 
ing to Eq. d2|) 



r = 9/uj{I). 



(157) 



It means that a minimum range, r, at which the tra- 
jectory has the coordinate and momentum equal to 
the given p and z, respectively, defines the desired 
angle variable, 9. 

In order to use Hamilton equations ( p9| ) and (J3C 
in practice we need a numerical procedure which al- 
lows one to express an arbitrary function F{p, z) as a 
function of I and 9, by replacing the arguments p and 
z with p{I,9) and z{I,9). In other words, we need 
a function $(1, 6) = F{p{I,9),z{I,9)). In principle, 
it can be found by a direct substitution. However, 
there exists a following more convenient procedure. 

Since both p{I, 9) and z{1, 9) are periodic in 9 with 
the period 2ir, the function <£>(/,#) is also periodic 
and it can be represented as a Fourier series 



$(J,0) = A n {I) cos n6 + ^B n (I) sin n8, 



n=0 



(158) 



where 



A n {I) 
B n {I) 
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d0$(I, 



cos n9 
sin n6 



(159) 



M n = 



1/2tt, n = 
1/tt, n > 



Using relation (157) we can express the above in- 
tegrals via the integrals over the ray trajectory com- 
puted with the ray code. This yields 



Mi) 

B n (I) 



2tv 



x / dr' F{p(r'),z(r')) 



M n u(I) 



cos ncu(I)r' 
sin nu(I)r' 



(160) 



If the function F depends not only on p and z, but 
on range r, as well, then the function $ and the co- 
efficients A n an d B n acquire an additional argument 
r. In Eq. (16C) this argument should be considered 
as a constant, i.e. there should be no integration over 
this argument. 

Note that functions p(r) and z{r) in Eq. ( |l60| ) 
which define a trajectory in the auxiliary waveguide 
may be quite different from range dependencies of 
coordinates and momenta for "real" ray trajectories 
satisfying Hamilton equations (|J) and (Q). 
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